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Abstract 

Here we present the results of a Large Eddy Simulation of 
a non-buoyant jet issuing from a circular orifice in a wall, 
and developing in neutral surroundings. The effects of 
the subgrid scales on the large eddies have been modeled 
with the dynamic large eddy simulation model applied 
to the fully 3D domain in spherical coordinates. The 
simulation captures the unsteady motions of the large- 
scales within the jet as well as the laminar motions in 
the entrainment region surrounding the jet. The com- 
puted time-averaged statistics (mean velocity, concentra- 
tion, and turbulence parameters) compare well with lab- 
oratory data without invoking an empirical entrainment 
coefficient as employed by line integral models. The use 
of the large eddy simulation technique allows examina- 
tion of unsteady and inhomogeneous features such as the 
evolution of eddies and the details of the entrainment 
process. 

Introduction 

Standard engineering design models for jets and plumes 
employ averaging techniques that mask the rich details 
evident in the real world. A snapshot of a plume billow- 
ing from a chimney or effluent discharging from an ocean 
outfall shows a highly irregular edge with regions of low 
concentration near the centre and pockets of high con- 
centrations fair downstream from the source. Superim- 
posed on this irregularity is a structure evidenced by the 
boils or eddies which seem to repeat at seemingly regular 
intervals. While experimental techniques are available 
to explore these unsteady phenomena, sampling is lim- 
ited. For example laser-doppler anemometry is a point 
method that misses the synoptic, while laser-induced 
fluorescence and particle-imaging velocimetry only pro- 
vide a slice of information with sampling limited to the 
frame rate of the video camera and the density of pixels 
in the charged couple device. An attractive alternative 
is to model the fluid motions at a sufficiently detailed 
scale so that the energy-containing eddies are resolved - 
this modeling approach is called Large Eddy Simulation 
(LES). We present here the results of such a study on 
a simple momentum jet discharging into a semi-infinite 
quiescent environment. 

The Model 

The equations to be solved are the filtered incompress- 
ible Navier-Stokes equations together with the continuity 
equation and a scalar transport equation representing a 
passive tracer (or pollutant). 

The non-dimensional governing parameters for this flow 
are: 

Reynolds number = Re = UqDo/u , 

Prandtl number = Pr = v/k , 


where Uo is the mean vertical velocity of the fluid leaving 
the source (used to non-dimensionalise all velocities), Do 
is the diameter of the orifice (used to non-dimensionalise 
lengths), u is the kinematic viscosity of the fluid, and 
k the scalar diffusivity. In the following we use T to 
represent the dimensionless scalar concentration (e.g. a 
tracer or pollutant). 

Large Eddy Simulation 

To create an LES computational model, the Navier- 
Stokes/continuity/ tracer equations are filtered with the 
filter cutoff chosen to remove the dissipative scales but 
retain the energy-containing scales. In the present study, 
the subgrid scales stresses and fluxes are modeled using 
the dynamic approach ([6], and [4]), which automatically 
determines, using different filter widths, the spatial dis- 
tribution of the magnitude of eddy viscosity as required 
by the subgrid-scale Smagorinsky model. This procedure 
obviates the need for any empirical determination of the 
Smagorinsky constant. 

A major advantage of this method is that it is much more 
likely to be successful for inhomogeneous flows, particu- 
larly in the present case where part of the domain is 
laminar. This issue was explored by Liu efc ai. [9] who 
compared the Smagorinsky and the dynamic models with 
particle image velocimetry laboratory measurements in 
the far field of a turbulent round jet. The Smagorinsky 
model was found to correlate poorly with the real turbu- 
lent stress, while the dynamic model yielded appropriate 
coefficients. The inadequacy of the Smagorinsky model 
was also demonstrated by Bastiaans et ai. [1] in an LES 
model of transient buoyant plumes in an enclosure. In 
this study, they found they could only get a match of 
the evolving plume by tuning the Smagorinsky constant. 
Such tuning reduces the generality of a model. 

Numerical scheme 

We use the LES code of Boersma, an adaptation of 
his own jet direct numerical simulation (DNS) code[3]. 
The code was modified by Basu [2] to compute buoy- 
ant plumes. A spherical polar coordinate system (R, 0, <f>\ 
along the radial, lateral and azimuthal directions respec- 
tively) is used here because, for the present flow with 
its conical mean growth, a spherical coordinate system 
allows for a well-balanced resolution of the flow field 
without excessive grid points. For presentation purposes, 
however, and for comparison to laboratory data, the re- 
sults are converted to the axisymmetric cylindrical coor- 
dinate system {z axial and r radial). The symbols U and 
V are used for velocities in the z and r directions respec- 
tively. Further details of the computational scheme can 
be found in [3], [4], and [6]. 

Boundary conditions and test parameters 


For the chosen staggered grid, only velocities normal to 
the boundary need be specified - i.e. six components. The 
equation for the tracer (T) is second order in each spatial 
direction, requiring six boundary conditions as well. 

For the bottom boundary a Dirichlet BC is used with log- 
law profile velocity at the orifice, to which white noise of 
2% rms is added. At the outflow boundary, or the top 
end of the current computational domain, we use the ad- 
vective boundary condition and at the lateral boundary 
we apply a stress- free condition ([3] and [7]). Periodic 
boundary conditions are used along the azimuthal 4> di- 
rection boundaries and lateral velocity (u#) is set to mean 
value at 6 = 0. 

The Reynolds number for the jet source is set at 3500 - 
within the range of experimental data for which the jet 
is turbulent at or near the nozzle[5]. 

The computations are carried out in a domain that ex- 
tends to 50 diameters downstream of the source. The do- 
main encompasses a conical volume of lateral angle tt/ 12, 
with a virtual origin that is 15 diameters upstream of the 
orifice. For the purpose of LES, this volume is discre- 
tised using a grid of size ( N r = 128, Ng = 40, N+ = 32). 
Time stepping is chosen to satisfy CFL conditions [3] 
which turns out to be typically 0.01 in dimensionless time 
units (Do/Uq). Spatially the grid is an order of mag- 
nitude smaller them the expected large scale (typically 
0.42, although there is an arbitrariness about the defi- 
nition of the jet edge) and an order of magnitude larger 
than the Kolmogorov length scale [10] as appropriate for 
LES modeling. The temporal sampling is much smaller 
than the Kolmogorov time scale. 

The computations are carried out for over 100,000 time- 
steps (a non-dimensional time of about t = 300). Aver- 
aging for statistical quantities was done towards the end 
of the simulations when we were convinced that the flow 
was statistically stationary. 

Model verification 

Confidence in the model can only come from comparison 
with laboratory data. There is a wealth of such data and 
we use a selection reviewed by List[8]. The data repre- 
sent air in air and water in water over a wide range of 
Reynolds number (but all in the turbulent regime). De- 
spite the range of fluids and laboratory conditions there 
is a surprising consistency between the results and hence 
we would expect our model to be able to match them 
well. 

To compare the model results, we have averaged first 
moment and second moment statistics over some 10,000 
samples. The particular sampling and averaging carried 
out yields an uncertainty of about 2% for any given sta- 
tistical estimate. 

Here we examine two statistics, the mean z-direction ve- 
locity ([/) and the tracer concentration (T). The pro- 
cedure for reducing the data was to fit Gaussian curves 
(using Matlab's routine nlinfit with emphasis on the in- 
ner radius) to the r-direction profiles. The fitted e-folding 
radius (r c ) and the fitted maximum value were used to 
scale the data. Of the many checks that we have made 
[11], we here present just two - we recall that the ex- 
perimental jet e-folding width grows linearly as 0.107z 
for U and 0.1 26z for T. (From the laboratory data the 
uncertainty for the growth rates is 3% for both U and T 
[ 8 |). 



Figure 1: Growth of jet width in terms of the e- 

folding width for z-direction velocity (circles) and tracer 
(crosses). 

The fitted slopes in the central region (beyond the ZFE 
and short of the apparent influence of the downstream 
boundary) are 0.105 for U and 0.126 for T (figure 1). 
The matches are remarkably close to the laboratory data 
and well within the data uncertainty. 

Other checks carried out [11] show that the centerline ve- 
locity and tracer decay as z~ l , and that turbulent fluc- 
tuations and turbulent momentum and scalar transport 
match laboratory data. In almost every case the match is 
well within data and model uncertainty. We also verified 
that most of the variance can be accounted for by the 
resolved portion of the model - as we would hope from 
an LES model. 

Unsteady features 

To examine some of the jet features we take sections both 
through the centerline (meridional) and across the jet (ra- 
dial). At the computer, it is instructive to roam through 
many such sections but even with just two we can see 
a great deal. Figure 2a is a meridional section between 
6 and 13 diameters from the source while figure 2b is 
a radial section at about 9 diameters from the source. 
We overlay the instantaneous tracer concentrations and 
velocity vectors. 

At the bottom centre of the meridional section (figure 
2a), between z = 6 and z = 8, we can see what at 
first sight appears to be a core region of almost source 
concentration with velocities matching the source veloc- 
ity. There is actually a peak concentration at z = 6.7, 
r = —0.2 that is at the center of a patch about to break 
from the core to form a more or less separate puff. Such 
puffs are shed repeatedly, and, in fact, the previously 
formed puff can be seen at z = 11. Given that the source 
is uniform and steady with the addition of a small level 
of noise it is remarkable how unsteady and irregular are 
the features. 

Within and about these large scale puffs are eddies of 
comparable scale. An anti-clockwise eddy with horizon- 
tal axis can be seen centered at z = 9, r = — 1. At its top 
this eddy transports jet fluid away from the core while at 
the bottom, ambient fluid is transported towards the jet 
center. The horizontal extent of the eddy and its three- 
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Figure 2: Instantaneous meridional (a) and radial (b) sections of tracer concentration and velocity vectors. The relative 
position of each section is indicated by a black line. 


dimensional character is evident from the radial section 
(Figure lb) at x = —1.5, y = -0.8. An example of an 
eddy with vertical axis can be seen in the radial section 
centered atx = 1.2, y = Q. As with the other eddy, 
this one is clearly advecting away core fluid in the out- 
ward portion and ambient fluid towards the center in the 
inward portion. Closer to the core, there appear to be 
several smaller scale eddies which further mix the jet fluid 
with ambient. Image sequences show that any jet fluid 
that strays too far from the center is readily re-entrained 
by the sink flow that exists outside the turbulent jet re- 
gion. 

Although we have been emphasizing the irregularity of 
the flow, there is an underlying apparent dominant size 
of the puffs and eddies, and a typical frequency at which 
the puffs are shed. In the following we look to quantifying 
these scales. 

Quantifying scale 

In an attempt to quantify the time scale we have ex- 
amined a time series of tracer concentration at a single 
point - on the centerline at 30 diameters from the source. 
This choice is arbitrary as any point should share the 
same temporal spectral characteristics. A raw spectrum 
is computed from a sequence of duration 460 time units 
at interval of 1.8 time units. Band averaging over 4 ad- 
jacent estimates reduces the uncertainty (figure 3). Al- 
though this still leaves broad confidence intervals, a clear 
peak is evident at 0.025 dimensionless frequency units 
(Uq/Dq). There may also be other peaks at about 0.07 
and 0.1 for which we offer three possible explanations - 
(1) they are just noise - note that we only use a 4- fold 
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Figure 3: Spectrum of tracer at z — 30. The dashed line 
is the raw spectrum, while the solid is band-averaged. 

band-average, (2) they represent real higher frequency 
motions - (3) they are spurious peaks associated with 
the non-sinusoidal shape of the puffs. We will have more 
to say about this with regard to the spatial spectrum. 

From the mean flow field one could invoke the Taylor 
frozen-field assumption to deduce a spatial scale, but an 
impediment to computing such a scale is the obvious in- 
homogeneity. However, we can make use of the scaling 
determined from the mean and variance to transform an 
axial sample of the tracer to a homogeneous equivalent. 







Figure 4: Spectrum of centerline scalar concentration in 
homogeneous transformed space 

Subtracting the mean centerline temperature and then 
dividing by the mean removes these two trends but still 
leaves a z- dependent wavelength. However, if we trans- 
form the z-axis by the antiderivative of the expected in- 
verse growth rate, the z-dependence is removed. In this 
case, the growth is expected to be linear, and hence the 
transformation is simply logarithmic. The spectra of 100 
such records were averaged to yield the spectrum shown 
in figure 4. 

The interpretation of the peak at 5 transformed wave 
number units is that the centerline eddy length is pro- 
portioned to distance from the source and is here equal 
to 0.2z. Curiously this measure of the large scale motion 
matches the mean jet width defined by the e-folding half 
- i.e. the energy containing eddies are the same order of 
magnitude as the jet itself. 

The secondary peaks at 6.2 and 9 sure reminiscent of those 
seen in the temporal spectrum (figure 3), except here, 
with the narrow uncertainty, we are sure they are not 
noise. Of the other two possible explanations (real higher 
wavenumber features or spurious peaks) we are convinced 
that the latter is the correct one. An examination of a 
centerline tracer record shows that the head of a puff is 
generally steeper than the tail. A spectrum computed 
from an artificially generated monochromatic wave with 
such an asymmetry recreates the spurious peaks seen 
here. 

Conclusion 

The LES of a simple jet does a remarkable job of sim- 
ulating laboratory experiment data without the need to 
invoke empiricism such as an entrainment coefficient. In- 
clusion of the large eddies in the model allows one to 
examine the detail of the entrainment process. The be- 
haviour of the jet is that of a series of puffs rather than a 
steady stream. The outward motion at the top of the puff 
together with the inward motion at the bottom acts to 
interchange jet fluid with ambient and smaller scale mo- 
tions. The axes of the eddies are not constrained, how- 
ever, to the lateral direction as eddies are also evident 
with axes at orthogonal orientation. A transformation of 
centerline tracer concentration yields an estimate of the 
large scales of the unsteady motions. For the case ex- 
amined here, this turns out to scale with distance from 


the source, and to be about the same size as the mean 
width of the jet. Whether this result is Reynolds number 
dependent remains to be seen. 
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